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Abstract 

We have used merger trees realizations to study the formation of dark matter haloes. 
The construction of merger-trees is based on three different pictures about the for- 
mation of structures in the Universe. These pictures include: the spherical collapse 
(SC), the ellipsoidal collapse (EC) and the non-radial collapse (NR). The reliabil- 
ity of merger-trees has been examined comparing their predictions related to the 
distribution of the number of progenitors, as well as the distribution of formation 
times, with the predictions of analytical relations. The comparison yields a very sat- 
isfactory agreement. Subsequently, the mass growth histories (MGH) of haloes have 
been studied and their formation scale factors have been derived. This derivation has 
been based on two different definitions that are: (a) the scale factor when the halo 
reaches half its present day mass and (b) the scale factor when the mass growth rate 
falls below some specific value. Formation scale factors follow approximately power 
laws of mass. It has also been shown that MGHs are in good agreement with models 
proposed in the literature that are based on the results of N-body simulations. The 
agreement is found to be excellent for small haloes but, at the early epochs of the 
formation of large haloes, MGHs seem to be steeper than those predicted by the 
models based on N-body simulations. This rapid growth of mass of heavy haloes is 
likely to be related to a steeper central density profile indicated by the results of 
some N-body simulations. 
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1 Introduction 



It is likely that structures in the Universe grow from small initially Gaussian 
density perturbations that progressively detach from the general expansion, 
reach a maximum radius and then collapse to form bound objects. Larger 
haloes are formed hierarchically by mergers between smaller ones. 
Two different kinds of methods are widely used for the study of the struc- 
ture formation. The first one is N-body simulations that are able to follow 
the evolution of a large number of particles under the influence of the mu- 
tual gravity from initial conditions to the present epoch. The second one is 
semi-analytical methods. Among these, Press-Schechter (PS) approach and its 
extensions (EPS) are of great interest, since they allow us to compute mass 
functions (Press & Schechter 1974; Bond et al 1991), to approximate merging 
histories (Lacey & Cole 1993, LC93 hereafter, Bower 1991, Sheth & Lemson 
1999b) and to estimate the spatial clustering of dark matter haloes (Mo & 
White 1996; Catelan et al 1998, Sheth & Lemson 1999a). 
In this paper, we present merger-trees based on Monte Carlo realizations. This 
approach can give significant information regarding the process of the forma- 
tion of haloes. We focus on the mass growth histories of haloes and we compare 
these with the predictions of N-body simulations. 

This paper is organized as follows: In Sect. 2 basic equations are summarized. 
In Sect. 3 the algorithm for the construction of merger-trees as well as tests 
regarding the reliability of this algorithm are presented. Mass-growth histo- 
ries are presented in Sect. 4, while the results are summarized and discussed 
in Sect. 5. 



2 Basic equations and merger-trees realizations 

In an expanding universe, a region collapses at time t, if its overdensity at 
that time exceeds some threshold. The linear extrapolation of this threshold 
up to the present time is called a barrier, B. A likely form of this barrier is as 
follows: 




(1) 
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In Eq.(l) a, f3 and 7 are constants, 5* = S*(t) = S 2 (t), where 5 c (t) is the linear 
extrapolation up to the present day of the initial overdensity of a spherically 
symmetric region, that collapsed at time t. Additionally, S = a 2 (M), where 
a 2 (M) is the present day mass dispersion on co moving scale containing mass 
M. S depends on the assumed power spectrum. The spherical collapse model 
(SC) has a barrier that does not depend on the mass (eg. LC93). For this model 
the values of the parameters are a = 1 and (3 = 0. The ellipsoidal collapse 
model (EC) (Sheth & Tormen 1999, ST99 hereafter) has a barrier that depends 
on the mass (moving barrier). The values of the parameters are a = 0.707, j3 = 
0.485, 7 = 0.615 and are adopted either from the dynamics of ellipsoidal 
collapse or from fits to the results of N-body simulations. Additionally, the non- 
radial (NR) model of Del Popolo & Gambera (1998) -that takes into account 
the tidal interaction with neighbors- corresponds to a — 0.707, (5 = 0.375 and 
7 = 0.585. 

Sheth & Tormen (2002) connected the form of the barrier with the form of the 
multiplicity function. They show that given a mass element -that is a part of 
a halo of mass Mq at time to- the probability that at earlier time t this mass 
element was a part of a smaller halo with mass M is given by the equation: 



ffa+za 4- \ac 1 \T(S,t/S ,t )\ 
f(S,t/S ,to)dS = -j= {AS)3/2 exp 



(AB) 2 
2AS 



dS (2) 



where AS = S-S and AB = B(S,t)-B(S ,t ) with S = S(M), S = S(M ). 
The function T is given by: 

T(S,t/S ,t ) = B(S,t) - B(S ,t ) + j2 [S °~ } S]n ^B(S,t). (3) 

n=l m ° b 



Setting Sq = 0, and B(So,to) = in Eq. 3, we can predict the unconditional 
mass probability f(S, t), that is the probability that a mass element is a part of 
a halo of mass M, at time t. The quantity Sf(S, t) is a function of the variable 
v alone, where v = 5 c (t)/a(M). Since 5 C and a evolve with time in the same 
way, the quantity Sf(S,t) is independent of time. Setting 2Sf(S,t) = vf{y), 
one obtains the so-called multiplicity function ][y). The multiplicity function 
is the distribution of first crossings of a barrier B{y) by independent uncor- 
rected Brownian random walks (Bond et al. 1991). That is why the shape of 
the barrier influences the form of the multiplicity function. The multiplicity 
function is related to the comoving number density of haloes of mass M at 
time t -N(M, t)- by the relation, 

„. , M 2 . .dmM 

"'<"> = ^^dhir (4) 
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that results from the excursion set approach (Bond et al. 1991). In Eq.4, Pb(t) 
is the density of the model of the Universe at time t. 

Using a barrier of the form of Eq.l in the unconditional mass probability, the 
following expression for f(u) is found: 

f(y) = v^A[l + Piav 2 )-^^)] exp (-0.5az/ 2 [l + P(av 2 )^} 2 ) (5) 



where 




V 

Fig. 1. Various multiplicity functions. The solid line corresponds to the spheri- 
cal collapse model (SC) of the PS approach. The dotted and the dashed lines are 
the predictions of the Ellipsoidal collapse (EC) and the Non-radial collapse (NR.) 
models, respectively. SC, EC and NR models are described by Eq.5 for different 
values of the parameters. Thick dashed and the dot-dashed lines correspond to the 
predictions of J01 and ST99 models, respectively. These models are described by 
equations 7 and 8, respectively. 
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In Fig.l we plot vf{y) versus v. The solid line results from Eq. 5 for the values 
of the parameters that correspond to the SC model. The dotted line is derived 
from the values of EC model while the dashed line corresponds to the NR 
model. For matter of completeness, we draw two more lines. The first line, the 
one with thick dashes- corresponds to the function of Jenkins et al. (2001, J01 
hereafter). This satisfies the following equation: 



In order to express the above relation as a function of z/, we substitute <t _1 (M) = 
v/5 c and we assume a constant value of <5 C , that of the Einstein-de Sitter Uni- 
verse, namely 5 C = 1.686. The above formula is valid for 0.5 < v < 4.8. The 
second additional line, dot-dashed one, is the function proposed by ST99 and 
is: 



where v' = v^fa and the values of the constants are: A = 0.322, p = 0.3 and 
a = 0.707. 

We have noted that the comparison of the above curves with the results of 
N-body simulations is in a very good agreement, except for the SC model. 
Such a comparison is presented by Yahagi et al. (2004). According to their 
results, the numerical multiplicity functions reside between the ST99 and J01 
multiplicity functions at v > 3 and are below the ST99 function at v < 1. Ad- 
ditionally, the numerical multiplicity functions have an apparent peak at v ~ 1 
-as those given by Eq. 5- instead of the plateau that is seen in the J01 function. 



3 The construction of merger-trees 

Let there be a number of haloes with the same present day mass M . The 
purpose of merger-trees realizations is to study the past of these haloes. This 
is done by finding the distribution of their progenitors (haloes that merged 
and formed the present day haloes) at previous times. A way to do this is by 
using Eq.2. First of all, one has to use a model for the Universe and a power 
spectrum. In what follows we have assumed a flat model for the Universe with 
present day density parameters Q m ,o = 0.3 and = A/3Hq = 0.7 where 
A is the cosmological constant and H is the present day value of Hubble's 
constant. We used the value Hq = 100 hKMs _1 Mpc _1 with h = 0.7 and a sys- 
tem of units with m unit = 1O 12 M /i _1 , r unit = l/i^Mpc and a gravitational 
constant G — 1. At this system of units H /H unit = 1.5276. 
Regarding the power spectrum- that defines the relation between S and M 



uf(u) = 0.315exp(- | 0.61 + Zn[<r -1 (M)] | 3 ' 8 ) 



(7) 



f(u) = A(l + is'~ 2p ) fi/^v' exp(-v' 2 /2) 



(8) 
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in Eq.2- we employed the AC DM formula proposed by Smith et al. (1998). 
The power spectrum is smoothed using the top-hat window function and is 
normalized for a 8 = a(R = 8/i _1 Mpc) = 1. 

For each one of the SC, EC and NR models, we studied four different cases. 
These four cases differ to the present day mass of the haloes under study and 
are denoted as SCI, SC2, SC3, SC4 for the SC model, EC1, EC2, EC3 and 
EC4 for the EC model and NR1, NR2, NR3 and NR4 for the NR model. The 
first case of each model (SCI, EC1 and NR1) corresponds to haloes with mass 
0.1 -measured in our system of units- while the second case (SC2, EC2 and 
NR2) corresponds to mass I. SC3, EC3, NR3 correspond to haloes with mass 
10 and SC4, EC4, NR4 correspond to masses 100. 

Pioneered works for the construction of merger-trees are those of LC93, Somerville 
& Kollat (1999), ST99 and van de Bosch (2002, vdB02). Many of their ideas 
are used for the construction of our algorithm, that is described below. 
First, let us define as N res , the number of realizations used. This number is 
the total number of present day haloes of given mass Mq. A mass cutoff M m i n 
is used, that is a lower limit for the mass of progenitors. No progenitors with 
mass less than M min exist in the merger history of a halo. M min is set to a 
fixed fraction of M . In what follows, M min = 0.05M (see also vdB02). Also, 
a min is the minimum value of the scale factor. We set a min = 0.1. Merging 
histories do not extend to values of a < a min . A discussion about the choice 
of a m i n is given in Section 4. 

Useful matrices which have been used are the following: M par (the masses of 
the parent haloes), M z (the masses that are left to be resolved) and M pro (the 
masses of the progenitors). The argument i takes the values from 1 to N max , 
where N max is the number of haloes that are going to be resolved. Initially, 
Nmax — Af res . The values of Aou, N res , M m i n , a m i n are used as input. We 
discuss the choice of the value of Au in Section 4. The scale factor a is set to 
its initial value, (that is the present day one), a = a = 1. Then, the matrices 
M par and M t take their values: M par (i) = M , M t (i) = M for % = 1, N res . It 
is convenient to set a counter for the number of time steps, so we set I step = 0. 
The rest of the procedure consists of the following steps: 
1st step: I step = I step + I. 

The equation 5 c (a p ) = Au + 5 c (a) is solved for a p , that is the new, (current) 
value of the scale factor. 
2nd step: for % — 1, N rnax 

3rd step A value of AS* is chosen from the desired distribution. 
4th step: The mass M p of the progenitor is found, solving for M p the equa- 
tion: AS = S(Mp) - S(M par {i)). If M p > Mi(i) then, we return to the 3rd 
step. Else, the halo with mass M p is a progenitor of the parent halo i . The 
mass to be resolved is now given by: Mi(i) = M t (i) — M p . If M p < M min then 
M p is too small and is not considered as a real progenitor, therefore we return 
to the 3rd step. Else, 

5th step: The number of progenitors is increased by one, (N pro = N pro + 1), 
and the new progenitor is stored to the list M pro (N pro ) = M p . 
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6th step: If the mass left to be resolved exceeds the minimum mass, (Mi(i) > 
M m in), we go back to the 3rd step. Else, we return to the 2nd step, where the 
next value of i is treated. 

7th step: We set the new number of haloes N max to be the number of pro- 
genitors, N max = N pro . We store the progenitors, M pro (j), j = l,N max . The 
masses of new parents are defined as the masses of progenitors by M par (j) = 
M pro (j),j = l,N max . The value of the scale factor is updated a = a p , as well 
as the mass left to be resolved, Mi(J) = M par (j),i = j,N max . If a < a min we 
stop. Else, we return to the first step. 

Output is the list of progenitors after a desired number of time steps. Re- 
garding the third step, we have the following remarks. For the SC case, the 
change of variables x = ^M-M^o) i eac l s to a Gaussian distribution with 

y/S(M)-S(M ) 

zero mean value and unit variance for the new variable x. So, we pick values 
x from the above Gaussian (that is the well-known and very fast procedure 
Press et al. 1990). The values of AS = xAoj are distributed according to the 
distribution described in Eq. 2. The distributions in the EC and NR cases 
cannot be expressed using a Gaussian, so a more general method is used. This 
method, described in Press et al. (1990) is as follows: 

Let's suppose that we would like to pick numbers from a distribution function 
G(x), where the variable x takes positive values in the interval [0, c]. First, 
we calculate the integral J c G(t)dt = A. We note that A is not necessarily 
equal to unity. Then, a number, let say y, between zero and A is picked from 
a uniform distribution. The integral equation: 

X 

J G(t)dt = y (9) 
o 



is solved for x. The resulting points (x, G(x)) are then uniformly distributed in 
the area between the graph of G, the x-axis and the lines x = and x — c. The 
values of x have the desired distribution. We have to note that this procedure 
is -as it is expected- more time demanding than the SC. This is mainly due to 
the numerical solution of the integral equation. However, it has the advantage 
of being general. 

In Eq.2, / becomes an one variable function for given values of So, to, an d 
t. Then, the above described general procedure is used in order to pick val- 
ues of S. For numerical reasons, it is convenient to use a new variable w = 

The distribution of the masses of progenitors is found using the following pro- 
cedure: Let N p be the number of progenitors of N res parent haloes at time t 
and let m m i n , m max be the minimum and the maximum mass of these pro- 
genitors, respectively. We divide the interval [m min , m max ] to a set of k max 
intervals of length Am. Then, the number of progenitors, N pro ^, in the in- 
terval [m m i n + kAm, m min + (k + 1) Am] for k — 0, k max — 1 is found. The 
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distribution of progenitors masses is then given by N m _ tree = A ^' k ■ This 
distribution has to be compared with the distribution that results from the 
analytical relations. Thus, multiplying Eq.2 by Mq/M, one finds the expected 
number of progenitors at t that lie in the range of masses M, M + dM. This 
number obeys the equation: 

N(M,t/M ,t )dM = ^f( S ,t/S ,t ) | £L | dM (10) 



which is known as the number-weighted probability for M. In Figs 2 and 3, 
the distributions of progenitors are plotted for all models studied, for two dif- 
ferent values of the scale factor. The results of Fig. 2 correspond to scale factor 
a = 0.375 and those of Fig.3 to a = 0.2. Details are given in the captions of the 




Fig. 2. The distribution of progenitors for various models used at scale factor 
a = 0.375 (redshift z = 1.67). All curves correspond to a halo with present day 
mass Mq = 0.1. The smooth dotted and the crooked dotted lines give log(N) 
and log(N rn ^tree) versus M/Mq for the SCI model, respectively. Smooth solid and 
crooked solid lines give log(N) and log(N m ^ t ree) for the NR1 model, respectively. 
In a similar way, the line with short thick dashes correspond to EC1 model. 
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Figs. Large differences between the predictions of spherical and non-spherical 
models are shown. These differences are very clear and they can be seen even 
if a relatively small sample of haloes is used (we used N res = 5000). On the 
other hand, the predictions for the two non-spherical models NR and EC are 
very close. 

The time evolution of halo's mass is known as its mass growth history. How- 
ever, a clear definition of the halo's mass at time t is needed. One can define 
as halo's mass at a step, the mass of its most massive progenitor at that step. 
We call this procedure the most massive progenitor approximation. A differ- 
ent approximation is the one followed by vdB02. We call this approximation 
the main progenitor approximation and it is as follows: At the first step, the 
most massive progenitor M mmp of a halo is found. This is the halo's mass at 




Fig. 3. The distribution of progenitors for various models studied at scale fac- 
tor a = 0.2 (redshift z = 4). All curves correspond to a halo with present day 
mass Mo = 10. The smooth dotted and the crooked dotted lines give log(N) 
and log(N m ^tree) versus M/Mq for the SC3 model, respectively. Smooth solid and 
crooked solid lines give log(N) and log(N m _ t ree) for the NR3 model, respectively. 
Lines with short thick dashes correspond to the results of EC3 model. The smooth 
line gives log(N) while the crooked one gives log{N m _tree)- 
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that time. At the next step, the most massive progenitor of the halo with 
mass M mmp is found and is considered as the new halo's mass. The procedure 
is repeated for the next steps. According to the definition of LC93, the scale 
factor at which the mass of the main progenitor equals half the present day 
mass of the halo is called the formation scale factor. We denote it by a/. A 
detailed description of the procedure follows: 

1st Step: I step — I step + 1- 

The equation 5 c (a p ) = Au + 5 c (a) is solved for a p that is the new (current) 
value of the scale factor. 
2nd step:for % = 1, N max 

The mass of the most massive progenitor of parent i is set equal to zero 
(M mmp (i) = 0). 

3rd. A value of AS from the desired distribution is chosen. 
4th step: A mass M p is found, solving for M p the equation: AS = S(M p ) — 
S(M par (i)). If M p > M t (i) then, we go back to the 3rd step. Else, the halo 
with mass M p is a progenitor. We set Mi(i) = Mi(i) — M p . If M p < M m i n then 
M p is too small and is not considered as a real progenitor and so we return to 
the 3rd step. Else, 

5th step: The most massive progenitor at this time step is found. 
M mmp (i) = max(M p , M mmp (i)). 

6th step: If the mass left to be resolved exceeds the minimum mass, (Mj(i) > 
M m in), we go back to the third step. Else, we check if the mass of the most 
massive progenitor is for the first time smaller than half the initial mass of the 
halo. If this is correct, then the formation scale factor of this halo is defined 
by linear interpolation: 

0.5M - M mmv [i) 
QforM = ap+ M par (i)-M mmp (i) {a - <+> 



Otherwise, we proceed with the next i. 

7th step: The list of the most massive progenitors is used, in order to find the 
mean mass at this time step. The time of the scale factor is updated a = a p , as 
well as the mass left to be resolved, Mi(j) = M par (j), j = 1, N res . If a < a m i n 
we stop, else we go to the first step. 

We have to note that the halo's mass at a step -derived by the above procedure- 
does not give the mass of its most massive progenitor at that step. For example, 
let a halo of mass Mo that has at the first step two progenitors with masses 
Mi and M 2 with Mi > M 2 . Then, the mass of the halo at the second step is 
defined as the mass of the most massive progenitor (Mi). This progenitor is 
not necessarily the most massive of all progenitors at that step, since it could 
be less massive than one of the progenitors of M 2 . The procedure of main 
progenitor described above has the advantage of being economic. It does not 
require the construction of a complete set of progenitors in order to find the 
most massive one at a step. 
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We tested both the most massive and the main progenitor approximations 
and we found small differences only for massive haloes. In order to be able to 
compare our results with those of other authors, we used the main progenitor 
approximation. 

A test for the reliability of merger-trees is the comparison of the distribution 
of formation scale factors with the one given by the analytical relations that 
follow. The probability that the mass of the progenitor, at time t, is larger 
than Mq/2 is given by: 

M 

P(t,M p > M /2)= J N(M,t/M ,t )dM (12) 

Mo/2 

The formation time is defined (LC93) as the epoch when the mass of the 
halo, as it grows by mergers with other haloes, crosses the half of its present 
day value. Then, P(t,M p > Mq/2) is the probability the halo with present 
day mass M had a progenitor heavier than M /2 at t, which is equivalent 
to the probability that the halo is formed earlier than t, P(< t,M ). Thus, 
P(t, M p > Mq/2) = P(< t, M ). Differentiating with respect to t, we obtain: 

^^2> = J ^L[N(M,t/M„,t )]dM (13) 

M /2 

which gives the distribution of formation times. 

The distribution of formation scale factors as it results from Eq.13 is plotted 
in Fig. 4. The left hand side snapshot corresponds to haloes with present day 




Fig. 4. The distribution of formation scale factors as it results from Eq.13 for the 
three models used in this paper. The dotted lines correspond to the SC model, the 
solid ones to the NR model and the dashed to the EC model. The left snapshot 
corresponds to haloes with present day mass 0.1 -in our system of units- while the 
right snapshot to haloes with present day mass 100. 
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Fig. 5. Distributions of formation scale factors derived either from analytical rela- 
tions (Eq.13) or from the analysis of merger-trees results (see Eq.ll). 

mass 0.1, while the right hand side one to haloes with present day mass 100. 
In every snapshot dotted, solid and dashed lines correspond to SC, NR and 
EC models, respectively. Notice that the distributions that correspond to NR 
and EC models are shifted to the left relative to the SC model. This means 
that structures formed earlier in the NR and EC models. As it can be seen in 
Fig. 4, the differences between spherical and non-spherical models are smaller 
for larger haloes. As it is known in the literature(e.g Yahagi et al, 2004 and 
references therein), the results of non-spherical models are in good agreement 
with the results of N-body simulations. Consequently, for massive haloes, the 
results of N-body simulations are close to the predictions of the spherical mod- 
els, as well. This disagrees with the results of vdB02 which predicts that the 
disagreement between spherical model and N-body simulations is larger for 
massive haloes. 

Fig. 5 consists of four snapshots. Bars show the distributions of formation 
scale factor as they result from Eq.ll. Wide bars are the predictions of SC 
model, while narrow filled bars are the results of EC or NR models. Solid lines 
are the predictions of Eq. 13 for the SC model, while dashed lines are the 
predictions of the same Eq. for EC or NR models. The two snapshots of the 
first row show the results of SCI and EC1 (left snapshot) and SCI and NR1 
(right snapshot) while the two snapshots of the second row show the results of 
SC4 and EC4 (left snapshots) and SC4 and NR4 (right snapshot). It is shown 
that the predictions of Eq. 13 are in very good agreement with the numerical 
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distributions. Additionally, it is clear that the differences between the distri- 
butions of various models, resulting from Eq. 13, are successfully reproduced 
by numerical distributions. It is also shown that differences between spherical 
and non-spherical models are more significant for haloes of small mass. 



4 Mass-growth histories 

We define as mass-growth history, MGH, of a halo with present-day mass Mo, 
the curve that shows the evolution of M(a) =< M(a) > /Mo where < M(a) > 
is the mean mass at scale factor a of haloes with present day mass M . Two 
input parameters are used in the method of construction merger-trees in the 
procedure described in the previous section. The first one is the 'time-step' 
Auj. For values of Auj < 0.3, it is shown (vdB02) that MGHs are not time-step 
dependent. According to that, we used a value Auj = 0.1. The second param- 
eter is a min . Our results are derived for a min = 0.1. In order to justify this 
choice, let us discuss the following points regarding the definition of the mean 
mass. Let N h (a) be the number of haloes at scale factor a that have masses 
greater than M m i n . Obviously, Nh(l) = N res . Let also M/ l) j(a) be the mass of 
the i th halo at scale factor a. We define by Mfc(a) the sum SMi,«( a )> where 
for k = 1 the sum is extended to all haloes that satisfy M/ l; j(a) > M m j n . For 
k = 2 the sum is extended to all haloes (obviously in that case, the masses of a 
number of haloes are equal to zero since, by the construction of merger-trees, 
their mass histories are not followed beyond M m i n ). Finally, for k = 3 the sum 
is over those haloes that satisfy Mh,i{a m i n ) > M m i n . Then, we can define the 
mean mass at a in three different ways: First M™ ean = Mi(a) / Nh(a) , second 
M mean = M 2 (a)/N res and third M 3 ~ = M 3 (a) / N h (a min ) . 
As a evolves from its present day value to smaller values, the number N h (a) 
becomes smaller and the three mean values defined above become different. 
It is clear that M™ ean overestimates the mean value, while M™ ean underesti- 
mates it. On the other hand, M™ ean requires a large number of N res so that 
the remaining sample of haloes N h (a min ) to be large enough. The choice of 
the value a min = 0.1 is justified by the fact that for a > a min the differences 
in the above three definitions are negligible. 

The following results have been derived using a number of 10000 haloes. We 
found that a number of haloes greater than 1000 is sufficient for producing 
smooth MGHs that are in exact agreement with the ones of our results. 
Fig. 6 consists of three snapshots. The first snapshot shows the results of the 
SC model. The other two snapshots correspond to the EC and to the NR 
model, respectively. Every snapshot contains four lines. From top to bottom 
these lines correspond to the cases 1, 2, 3 and 4, respectively. It is clear that -at 
present time- the slopes that correspond to smaller masses are smaller. This 
sequence of slopes is a sequence of formation times. Large haloes continue 
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iog(oc) iog(a') iog(a> 



Fig. 6. Mass growth histories for the three models studied. Lines - from top to 
bottom- correspond to cases 1, 2, 3 and 4 respectively. The growth rate of smaller 
haloes falls below some critical value earlier than that of larger ones. This means 
that smaller haloes are formed earlier 

to increase their masses with a significant rate, even at the present epoch, 
while the mass growth of small haloes is negligible. Thus, small haloes have 
been formed earlier. However, a formation scale factor can be denned by the 
condition that the growth rate of mass falls below some value. Bullock et al. 
(2001^1301) used as a characteristic scale factor, the one satisfying the relation 
diogM(a)/dlog(a) = //, where fi is a specific value, typically equal to 2. We 
used the value /i — 2 and we denote the solution of the equation a c . Addition- 
ally, we denote ay the solution of the Eq. M(a) = 1/2 and we plot the results 
in Fig. 7. This figure consists of two snapshots. Left snapshot shows log(a,f) 
versus log(M/M ) for the three models studied. From top to bottom, solid 
lines show the predictions of SC, NR and EC models, respectively. These lines 
are well fitted by linear approximations that are shown by dotted, dashed and 




log(M/M ) log(M/M ) 

Fig. 7. The relation between the 'formation' scale factor and mass. The left snapshot 
shows the relation between af and M. From top to bottom, the solid lines are the 
plots of log(a,f) versus log{M / Mq) for the SC, NR and EC models, respectively. 
The dotted, short-dashed and long-dashed lines are the linear fits of the above 
lines, respectively. From top to bottom in the right snapshot, the solid lines are the 
plots of log(a c ) versus log{M / Mq) for the SC, NR and EC models, respectively. The 
dotted line is the linear fit of the SC model, while the dashed line is a prediction of 
a c given by the toy-model of W02 (see text for more details) 
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long-dashed lines, respectively. Thus a/ and M are related by a power law 
a f (M) ~ M n . The values of n are n = 0.058, n = 0.074 and n = 0.085 for SC, 
NR and EC models, respectively. The solid lines in the right snapshot show 
log(a c ) versus log(M/M ). The dotted line is the linear approximation of the 
results of SC model. It can be seen that this approximation is very satisfac- 
tory. Linear approximations are not satisfactory for non-spherical models. The 
dashed line in this snapshot shows the predictions of the toy-model of Bullock 
et al. (2001, B01 hereafter) that calculates a c from the power spectrum using 
following procedure: Let a halo with mass M. The scale factor at the epoch 
of its collapse a c , is found by the solution of equation.: 

a[M*(a c )} = a(FM) (14) 

where, M* is the typical collapsing mass satisfying: <7[M*(a)] = 1.686D(1)/D(a), 
and F is a constant. D(a) is the growth factor predicted by linear theory 
(Peebles, 1980). The dashed line presented in the right snapshot of Fig. 7 is 
the prediction of the above procedure for F = 0.004 and it is a good approxi- 
mation of our results. 

A less accurate but very practical approximation of a c is given (see also Suto 
2003) by 

a c (M) « 0.23M - 13 (15) 

where, M is expressed in our system of units. A value of about 0.13 for the 
exponent of the above relation is confirmed by our results too, since the best 
linear fits of the solid lines in the right snapshot of the Fig. 7, are of the form 
a c (M) ~ M n with n = 0.1135, n = 0.12 and n = 0.1256 for the SC, NR and 
EC models, respectively. 

In a complementary paper of that of B01, Wechsler et al. (2002, W02 here- 
after), based on the results of their N-body simulations, proposed for the 
growth of mass of dark matter haloes the exponential model: 

M(a) = exp[-a c /i(l/a - 1)] (16) 

In Fig. 8, we compare our MGHs with the predictions of the model of W02. 
This figure consists of twelve snapshots. Every row has four snapshots that 
correspond to the four cases of masses (0.1, 1, 10 and 100) studied for every 
model. The predictions of SC model are presented in the first row. The sec- 
ond row shows the predictions of NR and the third one the predictions of EC 
model. The solid lines show our MGHs while the dashed lines are the predic- 
tions of the model of W02 given by Eq.16. The values of a c used are those 
shown by the solid lines of Fig. 7 that are predicted by the numerical solution 
of the equation dlog(M)/dlog(a) = 2. Thus, no other information (such as 



15 



the cosmology, power spectrum etc) is used but only the MGHs of haloes. All 
snapshots refer to M(a) > 0.01. 

The conclusion from Fig. 8 is clear. The predictions of Eq.16 are very good ap- 
proximations to the MGHs predicted by merger-trees realizations. The quality 
of fit depends on the mass of the halo and not on the model used. For example, 
the left column of Fig. 8 shows a very good fit for all haloes with mass 0.1. 
A more careful look at this figure shows that the fit becomes less satisfactory 
for the earlier stages of the evolution of large haloes. In fact, the predictions 
of merger-trees are steeper curves than those of the model of W02. This can 
be seen clearly at the 3thd and 4th column that show the predictions for 
haloes with masses 10 and 100. We note that the mass growth rate can be 
connected - under some assumptions- with the density profile of the halo (e.g. 
Nusser & Sheth 1999, Manrique et al. 2002, Hiotelis 2003). According to the 
'stable-clustering' model used by the above authors, haloes grow inside-out. 
This means that the accreted mass is deposited at an outer spherical shell 
without changing the inner density profile. We recall that the radial extent of 
a halo is defined by its virial radius R V i r) that is the radius that contains a 
mass with mean density A vir times the mean density of the Universe pb- A vir 
is, in general, a function of the scale factor a but in many applications it is 
set equal to a constant value of 200 (that is a value for A vir for an Einstein de 
Sitter Universe). The mass M vir contained inside R V i T at scale factor a satisfies 




Fig. 8. log(M/Mo) versus log(a) for all the models studied. The solid lines are the 
predictions of merger-trees realizations while the dashed ones are their fits by the 
model given in Eq.16 (see text for more details) 



16 



the two following equations: 



M vir (a) = 4tt / r 2 p(r)dr 



(17) 



o 



M vir (a) = -irA vir (a)p b (a)Rl ir (a) 



(18) 



If the mass history of the halo is known, (that is the evolution of its mass as 
a function of the scale factor), then the density profile can be calculated by 
differentiating Eq. 17 with respect to a, which gives: 



where dR V i r (a)/da is calculated by differentiating Eq.18. Finally Eqs. 19 and 
18 give, in a parametric form, the density profile. It is clear that the density 
profile calculated from the above approximation depends crucially on the mass 
growth rate. Hiotelis (2003) used a model based on the above assumptions for 
haloes that grow by accretion of matter. This means that the infalling matter 
is in a form of small haloes (relative to the mass of the growing halo) and 
the procedure can be approximated by a continuous infall. In that case, it was 
shown that large growth rate of mass (observed in heavy haloes) leads to inner 
density profiles that are steeper than those of less massive haloes (where the 
growth rate of mass is smaller). 

We recall that the density profile of dark matter haloes is a very difficult prob- 
lem contaminated by debates between the results of observations, the results 
of numerical simulations and the ones of analytical methods (e.g. Sand et. 
al 2004 and references therein). Especially the inner density profile seems to 
follow a law of the form p(r) oc r -7 but the value of 7 is not still known. Obser- 
vations estimate a value of 7 about 0.52 (Sand et. al 2004), while the results 
of N-body simulations report different values as 1 (Navarro et al. 1997), or 
between 1. and 1.5 (Moore et al. 1998). Additional results from N-body simu- 
lations, as well as from analytical studies, connect the value of 7 with the total 
mass of the halo (Ricotti 2002, Reed et al. 2005, Hiotelis 2003). In any case, 
it is difficult to explain the differences appearing in the results of different ap- 
proximations. For example, regarding N-body simulations, some of the results 
are effects of force resolution. In a typical N-body code, such as TREECODE 
(Hernquist 1987) the force acting on a particle is given by the sum of two 
components: the short-range force that is due to the nearest neighbours and 
the long-range forces calculated by an expansion of the gravitational potential 
of the entire system. As it can be shown, the value of the average stochastic 
force in the simulation is an order of magnitude greater than that obtained 
by the theory of stochastic forces. Consequently, small fluctuations induced 




1 



(a) da 




(19) 
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by the small-scale substructure are not 'seen'. This is the case of cold dark 
matter models in which the stochastic force generators are substructures at 
least three orders of magnitude smaller in size than the protoclusters in which 
they are embedded (e.g. clusters of galaxies). On the other hand, regarding 
the differences between the results of various approximations, we notice that 
some of these differences could be due to the presence of the baryonic mat- 
ter. Especially in regions, as the central regions of haloes, the large density of 
baryonic matter may affect the distribution of dark matter. 
Since in this paper we describe in detail technical problems related to the con- 
struction of reliable merger-trees, we have to note that the inside-out model for 
the formation of haloes described above, requires calculations for very small 
values of a, because the central regions are formed earlier. Therefore, in order 
to explore the center of the halo we have to go deep in the past of its history. 
Consequently, the problems relative to the calculation of mean mass that an- 
alyzed at the beginning of this section, become more difficult. The density 
profiles predicted by the above described approximation are under study. 



5 Discussion 

Analytical models based on the extended PS formalism provide a very useful 
tool for studying the merging histories of dark matter haloes. Comparisons 
with the results of high resolution numerical simulations or observations show 
that the predictions of EPS models are satisfactory (e.g. Lin et al. 2003, Ya- 
hagi et al. 2004, Cimatti et al. 2002, Fontana et. al 2004), although differences 
between various approximations always exist. These differences show the need 
for further improvement of various methods. 

In this paper we described in detail the numerical algorithm for the construc- 
tion of merger-trees. The construction is based on the mass-weighted proba- 
bility relation that is connected to the model barrier. The choice of the barrier 
is crucial for a good agreement between the results obtained by merger-trees 
and those of N-body simulations. As it is shown in Yahagi et al. (2004) mass 
functions resulting from the SC model are far from the results of N-body sim- 
ulations while those predicted by non-spherical models are in good agreement. 
Our results show that the distributions of formation times predicted by the 
non-spherical models are shifted to smaller values. Thus, the resulting forma- 
tion times are closer to the results of N-body simulations than the formation 
times predicted by the SC model. It should be noted that N-body simulations 
give, in general, smaller formation times than that of EPS models (e.g. Lin et 
al. 2003). It should also be noted that observations indicate a build-up of mas- 
sive early-type galaxies in the early Universe even faster than that expected 
from simulations (Cimatti, 2004). 

We have shown that a small sample of haloes (> 1000) is sufficient for the 
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prediction of smooth mass growth histories. These curves have a functional 
form that is in agreement with those proposed by the N-body results of other 
authors. In particular the set of relations: 



a c (M ) = 0.24M,' 



'0.12 



(20) 



and 



M(a)/M = exp[-2a c (M )(l/a - 1)] 



(21) 



describes very satisfactorily the results of non-spherical models (NR and EC). 
These relations give the mass M of dark matter haloes at scale factor a for 
the range of mass described in the text and for the particular cosmology used. 
M is the present day mass of the halo in units of lO 12 h _1 M . 
Finally, it should be noted that a large number of questions regarding the 
formation and evolution of galaxies remains open. For example, a first step 
should be the improvement of the model barrier so that the results of merger- 
trees fit better both the mass function (Yahagi et al. 2004) and the collapse 
scale factor (Lin et al. 2003). It should be noted also that structures appear 
to form earlier in N-body simulations and even earlier in real Universe. Ad- 
ditionally, predictions relative to the density profile are of interest. Some of 
these issues are currently under study. 
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